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HIGH ORDER RESOLUTION OF THE 
MAXWELL-FOKKER-PLANCK-LANDAU MODEL INTENDED FOR ICF 

APPLICATIONS 



ROLAND DUCLOUS, BRUNO DUBROCA, FRANCIS FILBET AND VLADIMIR TIKHONCHUK 



Abstract. A high order, deterministic direct numerical method is proposed for the nonrela- 
tivistic 2D X X 3D V Vlasov-Maxwell system, coupled with Fokker-Planck-Landau type operators. 
Such a system is devoted to the modelling of electronic transport and energy deposition in the 
general frame of Inertial Confinement Fusion applications. It describes the kinetics of plasma 
physics in the nonlocal thermodynamic equilibrium regime. Strong numerical constraints lead 
us to develop specific methods and approaches for validation, that might be used in other fields 
where couplings between equations, multiscale physics, and high dimensionality are involved. 
Parallelisation (MPI communication standard) and fast algorithms such as the multigrid method 
are employed, that make this direct approach be computationally affordable for simulations of 
hundreds of picoseconds, when dealing with configurations that present five dimensions in phase 
space. 



Keywords. High order numerical scheme, Fokker-Planck-Landau, NLTE regime, ICF, Magnetic 
field, Electronic transport, Energy deposition. 
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1. Introduction 




In the context of the interaction of intense, short laser pulses with solid targets |29j . Inertial 

Confinement Fusion (ICF) schemes [33] , the energy transport is an important issue. In this 
latter field of applications (ICF) , it determines the efficiency of plasma heating and the possibility 
to achieve the fusion conditions. The appropriate scales under consideration here are about one 
hundred of micrometers for the typical spatial sizes, and one hundred of picoseconds for the time 
scales. 

Several key features should be accounted for. First of all, in typical ICF configurations, a 
significant amount of energetic electrons have a large mean free path, exceeding the characteristic 
gradient length of the temperature and the density: the particles motion exhibits nonlocal features. 

A wide range of collisional regimes should be dealt with to describe the propagation and the 
deposit of energetic electrons from the underdense corona of the target to its dense and compressed 
core. 

The collisions are important even if the beam particles themselves are collisionless [3] : these 
particles, when propagating in a plasma, trigger a return current that neutralizes the inci- 
dent current. This return current is determined by collisions of thermal, background electrons. 
The structure of the generated electron distribution function is then often anisotropic, with a 
strongly intercorrelated two population structure. For nonrelativistic laser intensities, smaller 
than 10 18 Wcm -2 , a small angle description for collisions between the two populations is well- 
suited, leading to the classical Fokker-Planck-Landau collision model. The Coulomb potential 
involves a large amount of collisions with small energy exchanges between particles, so that the 
Landau form of the Fokker-Planck operator is required here. Such a configuration with two 
counterstreaming beams typically leads to the developement of microscopic instabilities that can 
modify strongly the beam propagation. We refer to the two-stream and filamentation instabili- 
ties, where the wave vector of the perturbation is respectively parallel and perpendicular to the 
incident beam [7J [8]. A self-consistent description of electromagnetic fields is then required to 
describe the plasma behaviour and associated instabilities. Furthermore in the process of plasma 
heating, strong magnetic fields are generated at intensity that can reach a MegaGauss scale and 
may affect the energy transport [BJ l2~Tl [23] . The sources of magnetic field generation include on 
the one hand the effects of the rotational part of the electronic pressure which is a cross gradient 
VT x Vn effect, and on the other hand the exponential growth of perturbations of anisotropic dis- 
tribution functions (Weibel instability). Some electromagnetic processes can be strongly coupled 
with nonlocal effects. 

The plasma model studied in this paper is based on the nonrelativistic Vlasov-Maxwell equa- 
tions, coupled with Fokker-Planck-Landau collision operators. It gathers the listed requirements 
at laser intensities which are relevant for ICF. At higher laser intensities, a relativistic treatment 
should be considered [3[ [36] , and collision operators with large energy exchanges are required if 
secondary fast electron production proves to be non-negligible, paticularly with dense plasmas. 
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There are several numerical methods that treat the Vlasov-Maxwell model together with 
Fokker-Planck-Landau type operators. Among them, the Particle-In-Cell (PIC) methods pro- 
vide satisfying results only in a limited range of collisional regimes. Moreover, they suffer from 
the "finite grid instability", that leads to numerical heating. Also the statistical noise and the 
low resolution of the electron distribution function by PIC solvers lead generally to an inaccurate 
treatment of collisions, particularly when dealing with low temperature and high density plasmas. 
Another approach consists in the expansion of the distribution function in Legendre polynomials, 
retaining the lowest order terms. However, with this approach, a strong anisotropy of the distri- 
bution function cannot be treated O [22j Q] . Both of these methods are well-suited in particular 
regimes but fail at modelling more complicated situations where a collisionless anisotropic fast 
electron population is coupled to a collision dominated thermal population. To overcome these 
difficulties, a spherical harmonic expansion has been proposed, that proves to be efficient [3]. 
Here we propose a different approach which consists in approximating the full model by a direct 
deterministic numerical method. It discretizes directly the initial set of equations and enables 
to preserve, at a discrete level, the physical invariants of the model (positivity of the distribu- 
tion function, total mass and energy, entropy decreasing behaviour, etc). Many deterministic 
schemes of this type have already been considered for homogeneous Fokker-Planck type opera- 
tors p~H El QUI G7] • The nonhomogeneous case, that includes the transport part (see [20] for a 
comparison between Eulerian Vlasov solvers), involves a large computational complexity that can 
only be reduced with fast algorithms. Multipole expansion [26] and multigrid [10] techniques, as 
well as fast spectral methods \27\ 119]. have been applied to the Landau equation. For computa- 
tional complexity constraints, very few results on the accuracy of these methods are known in the 
nonhomogeneous case |12[ 119] , particularly when the coupling with magnetic fields is considered. 

Our starting point for the transport part discretization is a second order finite volume scheme 
introduced in [12] . Its main feature is that it preserves exactly the discrete energy, if slope 
limiters are not active. We intoduce additional dissipation on these limiters in order to successfully 
address the two-stream instability test case. We will underline the important role of the limitation 
procedure for the accuracy, on the second order scheme. This scheme is compared in this test 
case with a fourth order MUSCL scheme [35], with a limitation ensuring the positivity of the 
distribution function [5j. A similar approach, with the introduction of a fourth order scheme for 
transport to avoid numerical heating, has already been proposed in the context of PIC solvers 
[32j . The discretization of the Maxwell equations is performed with a Crank-Nicholson method, 
allowing to have time steps of the order of the collision time. It is designed to preserve the discrete 
total electromagnetic energy, which is a very important numerical constraint when considering the 
coupling of Vlasov and Maxwell equations for applications aiming at capturing an accurate energy 
deposition. We use for the Landau operator a fast multigrid technique that proves to be accurate 
in a wide range of collisional regimes. Moreover, the use of domain decomposition techniques and 
distributed memory MPI standard on the space domain leads to affordable computational cost, 
allowing to treat time dependent 2D X x 3D V problems. As for the Lorentz electron-ion collision 
operator, we insist on discrete symmetry properties that are important when coupling to the 
Maxwell equations. 

Finally, we propose to validate the numerical method on several physical test cases. 

The paper is organized as follows. First, we present the model and its properties, then we discuss 
the numerical schemes for the transport part, their properties, and propose several numerical tests. 
Then the discretization for the collision operators is treated and we finally present physical test 
cases 1D X x 3D V and 2D X x 3D V that show the accuracy of the present algorithm. 

l 



2. Kinetic model 



(2) 



Two particle species are considered: ions which are supposed to be fixed (assuming an electron- 
ion mass ratio m e /mi <C 1), and electrons for which the evolution is described by a distribution 
function f e (t, x, v) where for the more general case (x, v) S ft x M 3 , with £2 C M 3 . The nonrela- 
tivistic Vlasov equation with Fokker-Planck-Landau (FPL) collision operator is given by 

(1) ~KT + V x • (v / e ) + — V v • ((E + V X B)/ e ) = C e , e (f e ,f e ) + C e>i (fe), 

ot m e 

where q e = —e is the charge of an electron and m e is the mass of an electron. On the one hand, 
electromagnetic fields (E, B) are given by the classical Maxwell system 

SE J 

— - c 2 V x xB = , 

Ot eo 

dB „ 

— + V X xE = 0, 
ot 

where eo represents the permittivity of vacuum, c is the speed of light. The electric current is 
given by 

J(i,x) = q e / f e (t, x, v) v dv. 
Moreover, Maxwell system's is supplemented by Gauss law's 

(3) V X -E=A V x -B = 0, 

where p is the charge density: 

p = q e (n e -n ) = q e I / e (i,x,v)dv - n 

and uq/Z is the initial ion density. 

On the other hand in (pQ), the right hand side represents collisions between particles, which 
only act on the velocity variable, so we drop the x variable. The operator C e)e (/ e ,/ e ) stands for 
the electron-electron collision operator whereas C ei i(f e ) is the electron- ion collision operator 

(4) C e MeJe) = / ^^ Vvf / $(-V-v')[f e (v')V v f e (-v)-f e (v)V v ,f e (V)]dV), 

whereas C e ^{f e ) is the electron- ion collision operator 

(5) C e>i (/ e ) = ^f 111 / V v • [d>(v)V v / e (v)] . 

o 7r eg m z e 

where In A is the Coulomb logarithm, which is supposed to be constant over the domain and ^(u) 
is an operator acting on the relative velocity u 

(a\ ^1 \ ll U H 2ld - U 

(6) *(u) = ||u||3 • 

As we assume ions to be fixed, the FPL operator can then be simplified for electron-ion collisions 
}12j . and reduced to the Lorentz approximation. We refer to [2] for a physical derivation. 

In this model, the Vlasov equation stands for the invariance of the distribution function along 
the particles trajectories affected by the electric and magnetic fields E and B. The Vlasov equation 
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representing the left-hand side in (fTJ) is written in a conservative form, but it can also be written 
in an equivalent non-conservative form, while Maxwell equations (J2])-(J3]) provide with a complete 
self-consistent description of electromagnetic fields. The coupling between both is performed via 
the Lorentz force term E + v x B in the Vlasov equation, and the current source terms in Maxwell 
equations. Furthermore, the FPL operator is used to describe elastic, binary collisions between 
charged particles, with the long-range Coulomb interaction potential. Classical but important 
properties of the system (pQ)-([3]) together with operators and ([5]), are briefly recalled. For 
detailed proofs, we refer to [12], [13] . 



2.1. Transport equation under electromagnetic fields. Let us neglect in this section the 
collision operators. The Vlasov-Maxwell system (pQ)-([3]) with a zero right-hand side is strictly 
equivalent to (JTJ)- (EJ) provided Gauss's laws ([3D are initially satisfied. This gives a compatibility 
condition at initial time. 

The mass and momentum are preserved with respect to time for the Vlasov-Maxwell system, i.e. 
system ([I])-© without collision operators 



^/ / e (i,x,v)( Mdxdv = 0, t>0. 



Moreover, conservation of energy can be proved for the Vlasov-Maxwell system by multiplying 
equation ([1]) by m e ||v|| 2 /2 and integrating it in the velocity space. It gives after an integration 
by parts 

f (e ||E(t,x)|| 2 + — ||B(t,x)|| 2 + f m e ||v|| 2 / e (i,x,v)dv 1 dx = 0, t > 0, 

with c 2 eo = 1 • The Vlasov-Maxwell system also conserves the kinetic entropy 

d f* 

-H{t) = - / / e (t,x,v)log(/ e (t,x,v))dxdv = 0, t>0. 

at JR 3 xR 3 



d 

Tt 



2.2. Collision operators. The FPL operator is used to describe binary elastic collisions between 
electrons. Its algebraic structure is similar to the Boltzmann operator, in that it satisfies the 
conservation of mass, momentum and energy 




= 0, t > 0. 



Moreover, the entropy is decreasing with respect to time 

^~(t) = fe(v,t) log(/ e (v,t))dv < 0, t > 0. 

dt dt J R 3 

The equilibrium states of the FPL operator, i.e. the set of distribution functions in the kernel of 
Ce,e(/e>/e)j are given by the Maxwellian distribution functions 



where n e is the density, u e is the mean velocity and T e is the temperature, defined as 

n e = / e (v)dv, 

U e = — / / e (v)vdv, 

n e Jr3 

Te = ^ [ / e (v)||v-U e || 2 dv. 

On the other hand, the operator ((Sj), modelling collisions between electrons and ions, is a Lorentz 
operator. It satisfies the conservation of mass and energy 



/ 



(',,(/, )(v)( i^.^na j dv •• (J. 



Moreover, the equilibrium states for this operator are given by the set of isotropic functions: 

Ker (C e ,i) = { f e € L 1 ((1 + ||v|| 2 )dv) , / e (v) = <j>(z), z = ||v - u e || 2 } . 
Finally, each convex function ip of / e is an entropy for C e j(/ e ), 



7 /> 

"77 / ^(/e)dv < 0, t >0. 



In addition to these properties, we present a symmetry property. This property may have some 
importance, in particular in presence of magnetic fields. In that case, any break of symmetry due 
to an inadequate discretization method could lead to generation of an artificial magnetic field, via 
the current source terms, when coupling with the Maxwell equations. 

Proposition 2.1. If f e has the following symmetry property with respect to the direction k at 
time to 

(7) / e (to,v) = / e (t ,V fc ), 

with components for 

v fc _ f +vj ifi^k, 
1 \ -v* if i = k. 
Then, this symmetry property is preserved with respect to time. 

3. Numerical scheme for transport 

We present a finite volume approximation for the Vlasov-Maxwell system ([I])-© without col- 
lision operators. Indeed, it is crucial to approximate accurately the transport part of the system 
to asses the collective behaviouiQ of the plasma, that occurs typically at a shorter scale than the 
collision processes. We introduce a uniform ID space discretization (xj +1 / 2 )j g /, I C N, of the 
interval (0, L\), in the direction denoted by index 1. The associated space variable is denoted by 
x\. We define the control volumes Cy = (a^-i/2) #1+1/2) x ( v j-i/2> v j+i/2)> the srze °f a control 
volume in one direction in space Ax and velocity Av. 



^By collective effects, we denote here the self-consistent interaction of electromagnetic fields and particles. Some 
collective effects are also considered in the collision processes, which make two particles interact via the Coulomb 
field. The self-consistent electromagnetic field then screens the long range Coulomb potential and removes the 
singularity in the Fokker-Plank-Landau operator. 
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The velocity variable v = t (vi,v 2 ,v 3 ) is discretized on the grid vj = j Av = t (vj 1 ,Vj 2 ,Vj 3 ) 
with j = *(ii, j 2 , j 3 ) G Z 3 . Moreover we note v j+1 / 2 = t {j 1 + 1/2, j 2 + 1/2, j 3 + 1/2) Av. Finally, 
the time discretization is defined as t n = nAt, with n £ N. 

Let /?j be an average approximation of the distribution function on the control volume Cy at 
time t n , that is 

$ - A~~m / /(*V,v)dxdv. 
' J Ax Av* Jc i j 

Moreover since the discretization is presented in a simple 1D X space geometry, the electromag- 
netic field has the follownig structure: E = t (Ei(t,xi)),E 2 (t,xi),0), B = *(0, 0, B 3 (t, x\)). Hence 
t (Ef i ,E2 i ) is an approximation of the electric field t (E\,E 2 ) whereas B 3i represents an approxi- 
mation of the magnetic field B 3 in the control volume (#i— 1/2, #1+1/2) a ^ time t n . 

3.1. Second order approximation of a one dimensional transport equation. For the sake 
of simplicity, we focus on the discretization of a ID transport equation; the extension to higher 
dimensions is straightforward on a grid, without requiring time splitting techniques between 
transport terms. In this section, the index 1 is dropped both on space and velocity directions, for 
this simple 1D X geometry. 

Let us consider the following equation for t > and x £ (0, L), 

(8) -dt +V dx- = ^ 

where the velocity v > is given. By symmetry it is possible to recover the case when v is 
negative. In the following we skip the velocity variable dependency of the distribution function. 
Using a time explicit Euler scheme and integrating the ID Vlasov equation on a control volume 

(^-1/2,^1+1/2), ^ Y ields 

(9) ft' = f? ~ |£ [-^+1/2 - 

where ^ + i/ 2 represents an approximation of the flux v f(t n , x i+ i/ 2 ) at the interface x i+1 / 2 . 

The next step consists in approximating the fluxes and to reconstruct the distribution func- 
tion. To this aim, we approximate the distribution function f(t n ,x) by fh{x) using a second 
order accurate approximation of the distribution function on the interval [^-1/2,^1+1/2), with a 
reconstruction technique by primitive [12] 

(x_ 
Ax 



(10) h(x) = f? + ef y —^(f? +l - fi 



We introduce the limiter 

r if (/f +1 - m <j? - fti) < 0, 



(ii) 



9 (\\ f°|| _ f n ) \ 

^(i, w fn _ fn Jt) ) if(/r +1 -/D<o, 

■li Ji+1 ) 
2f n \ 

min ( 1, — — — — I else, 
Ji+l ~ Ji 



and set ^ r ^_ 1 / 2 = v fh( x i+i/2)- This type of limiter introduces a particular treatment for extrema. 
At this price only (dissipation at extrema), we were able to recover correctly the two-stream in- 
stabililty test case, without oscillations destroying the salient features of the distribution function 
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structure. Another choice for the limitation consits in choosing the "Van Leer's one parameter 
family of the minmod limiters" |24j 



I fn fn\ I fn fn \ I fn fn s 

/ ln \ + • i / r Uj+1 Ji ) KJi+1 Ji—l) ,\Ji Ji-1, 

(12) ej = mmmod b — — , — - — , b- 



Ax 2Ax Ax 



where 



minmod(x, y, z) = max(0, min(x, y, z)) + min(0, max(x, y, z)) , (x, y, z) G 



and b is a parameter between 1 and 2. We will see on the two-stream instability test case the 
importance of the choice for limiters. 

Finally, this reconstruction ensures the conservation of the average and maximum principle on 



3.2. Fourth order transport scheme. We turn now to a higher order approximation (fourth 
order MUSCL TVD scheme) [35]. This scheme has also been considered in [5], in the frame of 
VFRoe schemes for the shallow water equations, where the authors proposed an additional lim- 
itation. Here we note that an optimized limitation procedure is possible in our case, breaking 
the similar treatment for both right and left increments, and taking advantage of the structure of 
the flux in the nonrelativistic Vlasov equation: the force term does not depend of the advection 
variable. 

For this MUSCL scheme, we only provide here with the algorithm for the implementation of 
this scheme and refer to [5], [35J for the derivation procedure of this scheme. 
The high order flux at the interface Xj+i/2, at time t n reads 



-pn _ -p I fn fn \ _ \ v '/i> ^ V > , 

^+1/2 - J- l/i,r> Ji+1,1) ~ | v jn_ if y < o 



This numerical flux involves the reconstructed states: /" r = /" + (Af)f and = /" + (Af) i , 
where (Af)f are the reconstruction increments. 

An intermediate state /*, defined by - (/™ r + /* + /^) = /" si introduced. It is shown in [5j 

that the introduction of this intermediate state preserves, provided the CFL condition is formally 
divided by three, the positivity of the distribution function. Following [35J and [5], the fourth 
order MUSCL reconstruction reads 
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Algorithm of reconstruction. 

Compute 

(A/)" = -i (2A*/ f _ 1/2 + A*/ m/2 ) , 
(A/)+ = i(A*/ i _ 1/2 + 2A*/ i+1/2 ) J 

where 

A*/i_i/ 2 = minmod(A*/i_ 1/2 ,4AVi+i/2)> 
A*/i+i/2 = mimnod(A*/ i+1/2 ,4A*/i-i/2) 

and 

1 3 - 

A /i+1/2 = i+1/2 — gA /i+1/2, 

A 3 /i+i/ 2 = A^ a _ 1/2 - 2A/f +1/2 + A/f +3/2 , 

with 

A A-i/2 = mmm od(A/ i _i/2,2A/ i+1 /2,2A/ i+ 3/ 2 ), 
i+1/2 = minmo d( A fi+i/2,2Afi+3/2,2Afi_ 1/2 ), 
A fi+3/2 = minmod(A /i+ 3 /2,2A/j_ 1/2 ,2Aj Vf-1/2), 
with the notation A/ i+1 / 2 = /i+i - /». 



Reminding that the minmod limiter is given by 



minmod(x, y) = < 
with (x, y) G M 3 . 



0, if xy < 0, 
x if |x| < |y|, 
y else, 



The limitation proposed in [5] is then applied. 
It allows to satisfy the positivity of the reconstructed states. 
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Algorithm for the limitation involving the intermediate state. 

Compute (Af) l ^ m,± such that 
f n + ( A /)f m - > 0, 

f n + (A/)f m ' + > 0, 

and 

/; = /r-(A/); im '--(A/): im ' + >o. 

This limitation reads: 



(a/); : 



lim, — 



max 



{(Af)r-fn), 



where 



(Af) l r + = 0max((A/)+, -f?) , 



( 1, if max ((A/)-, -/f ) + max ((A/)+ -/!») < , 



min ( 1 



£ n 

J i 



' max((A/)- ,-/" ) +max((A/)+ ,-f? ) 



otherwise. 



3.3. Application to the Vlasov-Maxwell system. We exactly follow the same idea to design 
a scheme for the full Vlasov equation in phase space (x, v) £ !] x I 3 , In addition, a centered 
formulation for the electromagnetic fields is chosen: 



(13) 



E" +1 / 2 = - (E" +1 + E") and B™ +1 / 2 = - (B 



n+l 



+ B n ) 



The discretization of the Maxwell equations (J2])-([3|) is performed via an implicit ^-scheme, with 
6 = 1/2, which corresponds to the Crank-Nicholson scheme and thus preserves the total discrete 
energy. This discretization is presented in a simple ID space geometry. The electric field E = 
t {E u E 2 , 0) and the magnetic field B = *(0, 0, B 3 ) are collocated data on the discrete grid. These 
fields are solution of the system 



(14) 



77171+1 



E 



l.l 



At 



Eli 



At 



+ c 



At 



jn 

J l,i 



R n+l/2 _ R n+l/2 ™ 
2 3,i+l "3,1-1 _ J 2,i 



2Ax 



rjn j-,71+1/2 
B 3,i f^,i+l 



n+l/2 
2,i— 1 



2Ax 



This scheme is well suited for the electrodynamics situations that are treated here in the test 
cases. 

n 



The approximation for the current in (|14p J" and J^ 1 has been chosen such as 



(15) Jli = £ Ay3 ^ /*j and J £ = E Ay3 ^ /y ■ 

Unfortunately, these expressions do not preserve the total energy when slopes limiters are active, 
but we will show that they have the important feature to reproduce the discrete two-stream 
dispersion relation. 

First, we remind discrete properties concerning positivity, mass and energy conservation [12] of 
the second order scheme (f9|)- (fT0j) coupled with (fT3|) - (fT5|) . considering now the magnetic component. 

Proposition 3.1. Let the initial datum (/fj)ijgz 3 be nonnegative and assume the following CFL 
type condition on the time step 



(16) 



At < C min ( Ax, Av) 



where C > is related to the maximum norm of the electric and magnetic fields and the upper 
bound of the velocity domain. 

Then the scheme <f^j)-< f 10\) coupled with A13\) - (T5\) . when extended to the infinite 3D X x 3D V 
geometry, gives a nonnegative approximation, preserves total mass and preserves total energy 
when slopes limiters are not active on the transport in the velocity directions 



- Ax- < 



e ||E?|| 2 + — ||B?|| 2 + m e 
Mo 



C° , n G N. 



In addition to these properties, we justifiy our choice for the numerical current thanks to a 
discrete dispersion relation on the two-stream instability. In the rest of the section, we drop the 
index 1 on the variables x\, v\, E\ and Ji, because the transport is considered 1D X x 1D V . 

Proposition 3.2. Consider the second order scheme IPjj-t fiOj) coupled with A13 \) -(T5 \) . when slope 
limiters are not active, to approximate the Vlasov- Ampere system 

dv 



(17) 



f df 


+ v 9£ 
dx 


dt 
< 


dE 


J 


> ~dt 





0. 



Then the definition Iil5\) for the current J defines a discrete dispersion relation that converges 
toward the continuous dispersion relation when Av, Ax and At tend to 0. 

Proof: The two-stream instability configuration can be fully analysed with the Vlasov-Ampere 
system (fTT|) extracted from equations dlj)-((3j) . The dispersion relation for a perturbation oc 
e i(kx-uit) Q £ an i n iti a j equilibrium state with ||/^|| -C ||/^||, then reads 



(18) 



1 + 



0/(0) 



eo f^e Jr w(u) — kv) dv 



dv = 0. 



Here the crucial point is the discretization on the velocity part of the phase space, so that we 
perform a semi-discrete analysis. In the frame of the discretization ([9])- (|10p coupled with (|13p -(|15p. 
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we consider the semi-discrete scheme approximating (|17p 



(19) 



( df 


+ v d£ 
dx 


at 


dE 


- Si 


~dt 


eo 



Av 



with 



fj+l/2 



fj+l + fj 



assuming the slope limiter is not active. Then we performe a discrete linearization around an 
equilibrium state 

where <C ||/(°)||. Using f (1) oc e i ( fca, - w *) in (HSJ>, it yields 



(20) 



- feuj) / + — E 

•> m e Av 



f 



(o) 



_ f(°) 

, ^(i) Jj+l/2 Jj-1/2 



-I ijJ 



These equations lead to the discrete dispersion relation 



(21) 



1 + 



AO) _ AO) 

Jj+l/2 Jj-1/2 

Av 



Av = 0. 



We recover the continuous dispersion relation (|18p when passing at the limit Av — > 0. Any other 
choice for the discrete current in (120p would introduce an additional error to the 0(Av 2 ) error in 
the relation dispersion (|21[) . For instance, choosing 

J = 5^Auu i / i+1/2 



would have lead to the analogous of (|21|) : 



(22) 



1 + 



E 



en m e ^— ' u> (oj — k v*) 



AO) _ AO) 

Jj+l/2 JJ-1/2 

Av 



Av = 0, 



which is a "shifted" dispersion relation, with a O(Av) accuracy, compared to the 0(Av 2 ) accuracy 
on relation (|21[) , □ 
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4. Validation of the transport schemes 



We first propose a validation stategy in the linear, collisionless regime, based on the work of 
Sartori and Coppa [30]. They performed a transient analysis, and obtain exact solutions of the 
periodic Vlasov-Poisson system, in the nonrelativistic and relativistic regime. 

Their approach, relying on Green kernels, is recalled in Appendix [A] in the nonrelativistic 
regime. A generalization of the 2D periodic relativistic Vlasov-Maxwell system, including mag- 
netic fields, will be presented in a forthcoming paper. Our objective is to capture kinetic effects in 
the linear regime, such as the Landau damping and the two-stream instability. A semi-analytical 
solutin is obtained, with a prescribed accuracy. Moreover, this method allows to explore wavenum- 
ber ranges where other approaches relying on dispersion relations fail. We recall that classical 
validations of kinetic solvers dedicated to plasma physics |12[ [25] are based on the calculation of 
the growth rates (instability), or decrease rates (damping) in the linear regime. Let us show the 
efficiency of the semi-analytical method on the two-stream instability test case. 



4.1. Scaling with plasma frequency. Scaling parameters can be introduced to obtain a di- 
mensionless form of the Vlasov-Maxwell- Fokker-Planck equations. The plasma frequency io pe , 
the Debye length Ad, the thermal velocity of electrons Vth, and the cyclotron frequency ui ce are 
defined as follows 



(23) 



n e z 



pe 



Vth 



eB 

me. 



€ m e v n e* 

These parameters enable us to define dimensionless parameters marked with tilde. 

• Dimensionless time, space and velocity, respectively: 

x v 

t = LO pe t, X = — , V = . 

• Dimensionless electric field, magnetic field and distribution function, respectively 



(24) 



eE 



eB 



(25) E = — , B 

This leads to the following dimensionless equations 

dfe 



0J, 



f _ f "th 
> Je — Je 



pe 



n 



(26) 



where (5 



dt 

dE 

~dt 

dB 

~dt 



+ V x • (v/ e ) — V v • ((E + V X B)/ e 



CeAfeJe) + uC e ,i(f e ), 



ff 2 



Vx x B 



n u, 



+ V x x E = 0, 



V x -E 



(1 - n), V x -B = 0, 
vth/c, v is the ratio between electron- ion collision frequency and plasma frequency 



Z no e 4 In A 



Z In A 

! An 



i 7T no A D OJ. 
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— with v e 

'pe 



Z no e In A 



The zero and first order moments of the distribution function are 

n(t,x) = / / e (i,x,v)dv, 



u 0, x ) = 77 x / /e(i,x, v) vdv. 
Moreover, in (j26[) the dimensionless collision operators are considered 

C e ,e(fe, fe) = V v • ( [ $(v - v') [/ e (v') V v / e (v) - / e (v)V v -/ e (V)] dv') , 

(27) { ^ Jr3 ' 

k Ce,i(fe) = V v • [*(v)Vy/ e (v)] , 

with $ given by ©. 

4.2. Test 1 : ID two-stream instability. The ICF physics involves a propagation of electron 
beams in plasma. The plasma response to the beam consists in a return current that goes opposite 
to the beam in order to preserve the quasineutrality. This leads to a very unstable configuration 
favorable to the excitation of plasma waves. We focus here on the instability with a perturbation 
wavevector parallel to the beam propagation direction, namely the two-stream intability. Of 
course, this stands as an academic test case but it is closely related to the physics of the ICF. 
Also it is a very demanding test for numerical schemes of transport, that have to be specially 
designed (see Proposition l3.2p . In particular, a discrete dispersion relation relative to that problem 
is developed to justify numerical choices for the second order scheme. For this scheme also, during 
the limitation procedure, an additional dissipation at extrema is introduced, compared to |12j . 
in order to preserve the solution from spurious oscillations. We will show the sensitivity of the 
scheme with respect to the chosen limiter, for this particular test case. Moreover, the fourth 
order scheme is introduced to reduce numerical heating, for simulations intended to deal with the 
two-stream instability. 

The (1D X x 1D V ) Vlasov- Ampere system (|17[) is approximated on a Cartesian grid. For this 
test case, we consider the scaling ()23j) - f)25j) . The initial distribution function and electric field are 

f(x,v) = -[{1 + A cos(kx))Mi, Vd (v) + (1- Aco8(kx))Mi- Vd (v)] , 



E°(x) 



0. 



where 



Mi, Vd (v) 



-^H 2 /2 



2vr 



is the Maxwellian distribution function centered around Vd- 

In order to compare the numerical heating associated with the second order and the fourth 
order scheme, we choose a strong perturbation amplitude A = 0.1. The perturbation wavelength 
is k = 2tt/L and the beam initial mean velocities are Vd = ±4, L = 25 being the size of the 
periodic space domain. We choose a truncation of the velocity space to be v max = 12 and time 
steps are chosen to be At = 1/200. 

The objectives of this numerical simulation are on the one hand to compare the second order 
finite volume scheme (specially designed to conserve exactly the discrete total energy, exept if the 
slope limiters are active) for different slope limiters and the fourth order MUSCL scheme. On 
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the other hand we want to explore the effect of a reduced number of grid points on the discrete 
invariants conservation. 
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Figure 1. Beams phase space (a) at initial time, (b) at 20 plasma periods (after saturation) 



In Figure [U two countersteaming beams that are initially well separated in the phase space (a) 
start to mix together. They finally create a complicated vortex structure, involving wave-particle 
interactions. This behaviour remains quantitatively the same whatever the transport scheme is 
(second or fourth order). However with a reduced number of grid points (smaller than 128 points 
in velocity), the second order (with limiter (fTTj) ) and fourth order schemes present a different 
behaviour for the total electric energy and total energy. 
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Semi-Analytical 



Dimensionless time 
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128 Points 
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Semi-analytical ~ 



20 25 30 35 

Dimensionless time 



(b) 



Figure 2. Evolution of the electrostatic energy for 32 2 ; 64 2 , 128 2 , 256 2 grid points, 
and the semi- analytical solution in the linear regime. Results are shown for (a) the 
second order with limiter (jlip and (b) fourth order transport scheme 
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Figure 3. Evolution of the electrostatic energy for 32 2 , 64 2 , 128 2 , 256 2 grid points, 
and the semi- analytical solution in the linear regime. Results are shown here for 
the second order scheme with the limiters (1121). with 6 = 2. 
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Dimensionless time 



20 30 40 

Dimensionless time 



(a) 



Figure 4. Comparison of the energy evolution for the second (with limiter (jlip ) 
and fourth order transport schemes. Results are shown (a) for 32 2 (6) 64 2 grid 
points 



For reduced grid resolutions, of 32 2 or 64 2 points, the fourth order scheme proves to be better 
than the second order one. For 32 2 points, plasma oscillations at the plasma frequency in the 
nonlinear phase are not reproduced with the second order scheme whereas they can be seen with 
the fourth order scheme (see Figure [2]). Moreover for this resolution, the transition from the linear 
phase to the nonlinear phase occurs earlier than it should for the second order scheme. 
As the grid resolution increases, the accuracy remains better for the fourth order scheme than for 
the second order one in the nonlinear phase (Figure [2]). The convergence toward curves with 128 2 
or 256 2 resolution grid is indeed better. We recall that quantities in Figure [2] and [3] are plotted 
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Figure 5. Comparison of the energy evolution for the second (with limiter (1121) . 
b = 2) and fourth order transport schemes. Results are shown (a) for 32 2 (6) 64 2 
grid points 



with a logarithmic scale, that smoothes out discrepancies between curves. In addition to these 
results, the respect of total discrete energy conservation proves to be better for the fourth order 
scheme than for the second order one at a reduced grid resolution, see Figure [4] and 

The use of limiters (112p for the second order scheme introduces accuracy improvements on 
the convergence behaviour and capture of plasma wave structure at reduce grid resolutions, see 
Figure [3j However, the energy dissipation remains quantitatively the same as the second order 
scheme with limiter (fTT|) , see Figures H] and [5l 

As this test case requires both a good preservation of invariants and accuracy when nonlinear 
phenomena occur, we might conclude that the fourth order scheme, with a resolution along each 
velocity direction greater than 32 cell, is well suited for our physical applications. The semi- 
analytical solution in the linear regime shown in Figure [21 using a Green function, brings some 
improvements compared to the classical validation in the linear regime, based on instabilities 
growth rates in the linear regime. In particular it discriminates precisely in time the linear and 
nonlinear phases. 

4.3. Test 2: ID X-mode plasma in a magnetic field. This test case stands as a validation 
in the linear regime for the coupling between Vlasov and Maxwell equations without collisions. 
A particular initial data is chosen (see the derivation in the appendix [B]) to trigger an X-mode 
plasma wave at a well-defined frequency u>. This type of wave presents a mixed polarization 
(longitudinal and transverse with respect to the magnetic field), that propagates in the plane P±, 
perpendicular to the magnetic field direction. 

The chosen frequency u is a solution of the dispersion relation (|7ip of the linearized Vlasov- 
Maxwell equations, introducing the equilibrium state /(°) (||v|| 2 ). The initial data are chosen 
such that f( ', E\, E2, and B3 only depend on u, ki = 2vr/Li and A; where f n , B 3 , E ± and 

Ei are the reconstructed (in the appendix [B]) Fourier transforms of the distribution function and 
electromagnetic fields. The magnetic field B' ' is the nonperturbed magnitude of the magnetic 
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field, L\ is the length of the space domain, A is the perturbation amplitude. The initial data can 
then be constructed with the help of truncated Fourier series 



/<°>(si,v) = / (0) (||vf) 



+ 



2 

£ 

n=-2 



/n(v±; 



J,k\ xi + i nip 



xi £ (0,Li), v G 



E 1 (t,x 1 ) = E ie - Mk ^\ Xl e (0,Li) , 
E 2 (t,x l ) = E 2 e- Mk ^\ xi€(0,L!) , 



■io;t+ifcia;i 



xi e (o,Li 



We define ip as the angle in the cylindrical coordinates for the velocity, defined with respect to 
the direction of the magnetic field (See appendix [B]) , 

The normalisations are defined by relations (|23p - (|25p . We choose = 2 and a rather strong 
amplitude perturbation A = 0.1 with periodic boundary conditions on the space domain. Also we 
have set (3 = Vth/c = 0.05. The dispersion relation have been solved for these parameters. One of 
the solution uj ~ 5.1432 is injected in the initial data set. 

We considered 126 points along the ID space direction, and 64 points along each velocity 
direction v = t {v\, x>2, v^). The dimension of the space domain is L\ = 25 whereas the truncation 
of the velocity space occurs at v max = 7 for each velocity direction. Furthermore, the time step 
is At = 1/200. 

The Fourier spectrum in Figure [6] exhibits a well defined frequency / = 1/T ~ 1.6375 (corre- 
sponding to a period T) for the total magnetic energy, that corresponds to a frequency f /2 for 

the magnetic field oscillations. We finally find uj = 2tt (t ) ~ 5.1443 from the numerical solution, 



to be compared with the analytical results 5.1432. This proves a good accuracy of the numerical 
results, while the distribution function is greatly affected by the magnetic field. As an illustration, 
we show in Figure [7] how the magnetic field makes the distribution function rotate in the velocity 
space perpendicular to the magnetic field axis. 



5. Approximation of the collision operators 

In the following, the presentation is restricted to the space homogeneous equation, for the sake 
of simplicity, 

, /(0, v) = /(°)(v), 
where C e>e (f, f) and C e>i (f) are given by §27$. 

5.1. Discretization of the Lorentz operator. We consider /j an approximation of the distri- 
bution function /(vs) and introduce the operator D, which denotes a discrete form of the usual 
gradient operator V v whereas D* represents its formal adjoint, which represents an approxima- 
tion of — V v -. Therefore, for any test sequence (^j)j e z 3 > we se t (D^jez 3 as a sequence of vectors 
ofM 3 
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Figure 6. Discrete Fourier spectrum in frequency of the discrete analogous of the 



total dimensionless magnetic energy 



Ll IIS*" 2 



-dx\. 



where D s is an approximation of the partial derivative -Mp with s £ {1,2,3}. In order to preserve 
the property of decreasing entropy at the discrete level, we use the log weak formulation of the 
Lorentz operator pi 



C e ,<(/)(v) V(v) dv = - / $(v) /(v) V v log(/(v)) • V v V(v) dv, 

where <3? is given by © and V is a smooth test function. Then, using the notations previously 
introduced, the discrete operator C^(f) is given by 



(28) 



a A r(/)(vj) 



-D' 



r.P 



5(v j )/ j D(log(/ j )) 



where S'(Vj) is the following matrix 



5(vj) 



|vj || 2 Id — vj (%> Vj. 
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Figure 7. Projection on the v\ — V2 velocity domain for the distribution function 
is shown at initial time t n = and at t n = 18.72, for a particular point of the 
space domain, x\ = 23.0114, V3 = 0. 



Now, vj has to satisfy the discrete conservation of energy 

(29) ^idlvjll 2 ) = ^2(||vj|| 2 ) = ^3(||vj|| 2 ) 

Vji v i2 v i3 

Then, we consider the 8 uncentered operators D e , with the formalism: 

D £ = \D?,D?,D«), 

with e = *(ei,£2,e3), and e, G {+1,-1} for i £ {1,2,3}. More precisely, the operator D €i is the 
forward uncentered discrete operator if = +1 and the backward uncentered discrete operator 
if £,: = -!: 



(30) D e * 




h- 


f-ei 




32- 






33- 







This 8 operators respectively match to 8 expressions of v| , following ([29 



This choice has been made to avoid the use of the centered discrete operator that conserves 
non physical quantities. On the other hand, the uncentered operators, taken separately, intro- 
duce some artificial unsymmetry in the distribution function leading to a loss of accuracy when 
coupling to Maxwell equations. To overcome these difficulties, following the idea of [9], we intro- 
duce a symmetrization of the discrete operator based on the averaging over the eight uncentered 
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discretizations: 



CUf) 



D* 



^(vn/jD^iogC/j)) 



This final expression will introduce an additional discrete symmetry property compared to the 
operator presented in [12] , 

We now present the discrete properties for the electron-ion collision operator. We have the 
classical properties: mass and energy preservation, an entropy decreasing behaviour, the positivity 
preservation of the distribution function in a finite time sequence. The proofs are not detailed 
here but can be deduced easily from those presented in [12j . The difference stands in the fact 
that we obtain the operator as an average over the full set of the uncentered operators (instead 
of an average over two operators). This modification allows to get a discrete analogous of the 
symmetry property presented in Proposition l2.lt 

Proposition 5.1. Under the condition \29\) on Vj, the discretization $31\) to the Lorentz operator 
)[5j) satisfies the following properties, 

• it preserves mass and energy, 

• it decreases discrete entropy 



H(t) = Av 3 £/j(t) log(/j(t)), 



there exists a time- sequence At n such that the scheme 

= A" • A/r, A /i/;;v 



en 



ill 

+ 00. 



defines a positive solution at any time i.e. ^2 n t r , 
Furthermore, if /j is symmetric with respect to in the direction jk a t time t n , then this property 



is preserved at time t n+1 , 
(31) 



E C^(/)(vjK fc A* 
jez 3 



Proof: We prove the last property and rewrite the operator (|31f) in a different manner, assuming 
we have a symmetry along the velocity direction Vj k 



1 



(32) (vj) = g2j^(/)(vj) 

where the notation e^'^ refers to 
(33) 



±,(fc) 



±,(fe) 



±1 

e,; 



if i = k, 



,(fe) 



vj) + C, 



(AO/ v 
.., ( V j) 



22 



We are interested in the cancellation of the operator C^- J (/)(vj)uj fc . This is equivalent to the 
cancellation of 

Q« := E(^r(vj) + ^' (fc) (vj 



E 



■ 5( ^.« )D ^w log(/ . 



D v h 



3 h v j 



-• w )D e - ,W log(/j) 



• D u 



Then, since D c 



h(*0 



% = D e = e k , it yields 



Q 



(k) 



l 



- E 



1 



,3 H V JI 



.(*>) 



E ct/j E 



je: 



?3 il v J 



E 



,(fc) 



^' (fc) (log(/j)) 



x-^ +,(fc) 

E^: £ e < (log/j) 



Jz 



■(*) 



E 



^ T ' (fc) (log(/j)) 



^ (log/. 



Then using definition (|33p and the symmetry of /j 1 with respect to in the velocity direction Vj k , 
we obtain = 0. Then multiplying (|32|) by and integrating in the full velocity space gives 
the relation (|3ip . This relation implies that f^ +1 is symmetric with respect to in the direction 



5.2. Discrete Landau operator. We consider the discretization of the FPL operator (JH) on the 
whole 3D velocity space. It is based on the entropy conservative discretization introduced in |14j . 
where a discrete weak log form of the FPL operator is used. This discretization yields: 



(34) 



r am 

dt 



p(t) = Av 3 E /j(*)/m(t)*(v j -v m )(D(log(/(t)) j -D(log/(t)) I 



where D stands for a downwind or upwind finite discrete operator approximating the usual gra- 
dient operator V v . This uncentered approximation ensures that the only equilibrium states are 
the discrete Maxwellian. The use of centered discrete operators would have lead to non physical 
conserved quantities. The discretization of the FPL operator is then obtained as the average over 
uncentered operators, but here for a different reason as in the previous section, on the electron-ion 
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collision operator discretization. In [10J, the scheme is rewritten as the sum of two terms: a sec- 
ond order approximation and an artificial viscosity term in Av 2 which kills spurious oscillations. 
However the computational cost of a direct approximation of (|34p remained too high. Therefore, a 
multigrid technique has been used. We refer to |10j and [11] for the details of the implementation 
on the FPL operator. Nevertheless, these latter approaches introduce a new approximation than 
can affect accuracy. Based on |27j . Crouseilles and Filbet proposed another approach and noticed 
that the discrete FPL operator (|34p in the Fourier space can be written as a discrete convolution, 
which directly gives a fast algorithm. Here we adopt the multigrid method, detailed in [TO], that 
has a complexity of order 0(n^ logra^). 

This discrete approximation preserves positivity, mass, momentum, energy, and decreases the 
entropy. Moreover the discrete equilibrium states are the discrete Maxwellian. 



6. Numerical results 



6.1. Scaling with collision frequency. For the analysis of collisional processes, a new scaling 
is introduced here, that allows time steps to be of the order of the electron-ion collision time. In 
order to account for transport phenomena occuring at the collision time scale, several parameters 
are required: the electron-ion collision frequency i/ e j, the associated mean free path A e j, the 
thermal velocity vth, and the cylotron frequency uj ce 



(35) 



Z no e 4 In A 
87re§rof i& 



Vth 



kbTq 



eB 



These parameters enable us to define the dimensionless parameters with tilde. 

• Dimensionless time, space and velocity, respectively 

(36) t = u eji t, x = x/X e>i , v = v/v th . 

• Dimensionless electric field, magnetic field, and distribution function, respectively 



eE 



eB 



(37) E = , B 

m e vth,v e ,i m e v e ^ 

This leads to the following dimensionless equations 



^ce f _ f V th 
i Je — Je 



( df 

-j£ + V x • (v/ e ) - V v • ((E + v x B)/ e 



(38) 



^Ce,e(/ej f e ) + C ei j(/ e ), 



<9E 
~dt 



1 



1 



^V x x B = - 2 nu, 



dB 

at 



V x -E 



V X B 



a" 
0. 



n) 



where a = u ei i/oj pe and (3 = vth/c. The collision terms C e;B (f e , f e ) and C ei i(f e ) are given in ([2?]) , 
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6.2. ID temperature gradient test case. In the context of laser produced plasma, the heat 
conduction is the leading mecanism of energy transport between the laser energy absoption zone 
and the target ablation zone. 

In such a system, the parameters of importance for the heat flux are 

• The effective electron collision mean free path A e . 

• The electron temperature gradient length At- 

• The magnetic field B and its orientation with respect to VT. 

These parameters enable to distinguish different regimes of transport, according to the Knudsen 
and the Hall parameters. 

The Knudsen number K n is a mesure of the thermodynamical non-equilibrium of the system 
(39) K n = ^. 

AT 

A regime characterized by K n — > refers to an hydrodynamical descripion, whereas a regime 
characterized by K n > 1 refers to a kinetic description, where the nonlocal phenomena appear. 
The parameters for ICF imply K n > 0.1, while the classical, local approach fails at K n > 0.01. 
This premature failure of the classical diffusion approach in plasma is explained by a specific 
dependence of the electron mean free path on their energy. In our applications the energy is 
transported by the fastest electrons, which have a much longer mean free path. 
The Hall parameter x = ^cT quantifies the relative importance of magnetic and collisional effects. 
oo c = eB/m e is the electron cyclotron frequency and r the mean electron- ion collision time 

^2 /^-^ 3 /2 



3 167r 2 e^/m^Te 



(40) r = - 

4 ^ ni Z 2 e 4 lnA 

For this test case, a simple gradient temperature configuration is shown in figure ([8]), modelling 
the following situation: through a layer of homogeneous plasma, a laser deposits its energy on the 
hot temperature side and the absorbed energy is transported with electrons to the cold tempera- 
ture side. 

Let us define the average over velocity of a function A(v) 

(41) < A >= — f Afdv , 

n e Jr3 

where n e (t, x) = / /(t, x, v)dv is the density of electrons. 

Jm. 3 

Following [61 [16] , we introduce the macroscopic quantities 

' j = q e n e (v) , 



(42) 



q = ^m e n e ((v • v)v) , 



R= / m e vC ei i(/ e )dv, 
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Figure 8. Initial configuration for the temperature gradient test case: a tem- 
perature profile is considered between to two domains of plasma with particles at 
thermodynamical equilibrium. Zero current boundary conditions enable to main- 
tain mass conservation. A heat flux is generated wherever there is a nonzero tem- 
perature gradient, as well as boundary layers on the heat flux, temperature, and 
electromagnetic profiles. 



(43) 



( 1 

p = n e T e = -m e n e ((v- < v >) • (v- < v >)) 



II = -m e n e ((v— < v >) <g) (v— < v >)) — pi, 



qioc = ^m e n e ([(v- < v >) • (v- < v >)] (v- < v >)) 



There, j is the electric current, q the total heat flow, R the friction force accounting for the 
transfer of momentum from ions to electrons in collisions, T e is the temperature, p is the scalar 
intrinsic pressure, II is the stress tensor, qi oc is the intrinsic heat flow and I the unit diagonal 
tensor. 

Quantities p, II and qi oc ar e defined in the local reference frame of the electrons, whereas j, q 
and R are defined relative to the ion center of mass frame. Ions are supposed to be at rest. We 
have the relation 
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(44) qi oc = q+j • (-pi + II)/(n e e) + j(-m e n e < v > 2 )/(n e e). 

The validation of our Fokker-Planck solver in the domain close to the hydrodynamical regime 
(local regime) requires knowledge of transport coefficients. Following the formalism of Braginskii 
[6] for the transport relations, the transport coefficients in the hydrodynamical regime have been 
calculated by Epperlein in [16]. These coefficients a ep , /3 e p> ^ep> are the electrical resistivity, 
thermoelectric and thermal conductivity tensors, respectively. From these quantities, we are able 
to compare the heat flux and electric field issued from the Fokker-Planck solver to those calculated 
analytically in [16], in the local regime. 

The classical derivation procedure to obtain the transport coefficients involves the linearization 
of the Fokker-Planck-Landau equation, assuming the plasma to be close to the thermal equilib- 
rium. The distribution function is approximated using a truncated Cartesian tensor expansion 
/(t,x, v) = /^(||v|| 2 ) + - 2 • fW(t,x, v). Following [IB], II and m e n e < v > 2 are neglected. 

Then considering appropriate velocity moments of f^ 1 ), electric fields and heat fluxes are expressed 
as a function of thermodynamical variables. The coefficients of proportionality, in the obtained 
relations, are defined as the transport coefficients. 

Several notations can be used, depending on the chosen thermodynamical variables. Adopting 
the Braginskii notations, we obtain 



(45) 
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We want to compare of the results of the solver with the analytical electric fields and heat fluxes in 
the local regime. For that purpose, we use the values of coefficients, for Z = 1, that are tabulated 
in [16]. As for the components of these tensors, we make use of the standard notations ||, _L, and 
A. Directions denoted with || and _L are respectively parallel and perpendicular to the magnetic 
field. Consequently, the parallel and perpendicular components of a vector u are respectively 
uu = b(u • b) and u±_ = b x (b x u), where b is the unit vector in the direction of the magnetic 
field. The direction defined by the third direction in a direct orthogonal frame is denoted by A. 
In the system ()45|) . the relation between any transport coefficient tensor ip and vector u is defined 
by 



(46) 



ip ■ u = y?iib(b • u) + C/2j_b x (u x b) ± (f/^b x u 



where the negative sign applies only in the case ip 
dimensionless form 



a ep . These coefficients can be expressed in 
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The dimensionless transport coefficients a^, /3£ , n^ p are functions of Z and the Hall parameter 
X = uj c t only. 

The heat flux and the electric field in (|45[) can then be rewritten in terms of dimensionless 
quantities, for the particular ID geometry of our temperature gradient configuration. In that 
case, the normalizations using a collision frequency ()35|) - (|37[) are used. 



(48) { 



q\ 
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-T e n- l j x - xT e B^V xl T e K c ep ^ - T e [P^ ± jt - 0^ A j 2 ) , 
5 

-T e U~ h ~ XT e B^ l V Xl T e Kl vA - T e (/?ep,±j2 + /?ep,Ail) . 



Ei=n e 1 j2B 3 -n e 1 V Xl p-V Xl T e (3^ p± + n e 1 B 3 x 1 (at P:± ji + al pA j 2 ), 
E 2 = -n e ~ x jxB 3 - V Xl T e Pep, A + ra e _1 ^3X _1 (aep,±J2 - "ep.Ail)- 
The Hall parameter x is expressed in terms of the dimensionless quantities B 3 and T e : 
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6.2.1. Electron transport in the local regime. In order to validate the numerical scheme in the 
local regime, we compare the heat flux Qfp and electric field Efp computed from the numerical 
solution, with those analytically (denoted by Qbr and Ebr) computed from the system 
The transport coefficients a ep , f3 ep , K ep have been tabulated in [16J. 



In this test case the source term can be considered stiff; the discretization of the collision 
operator is then of crucial importance and its accuracy can be tested. Moreover we provide, in 
this local regime, with validation results for a wide range of Hall parameters corresponding to 
ICF applications. 

The initial temperature gradient T e {x\) has the form of a step 



(50) T e (si) 




if x\ > xf , 
else , 



where T R and are third order polynomials in x\ — X™, x\ standing for the space coordinate 
and x™ for the mid-point of the ID domain. The coefficients of these polynomials are chosen such 
as they verify the following conditions at x™ 



(51) 



dx x y w d Xl v 1 ; (xf -x[)/X' 
T e L (x?)=T e R (x?) = ^±^, 
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and at the boundaries 




(52) 
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dT e R 




where Tl (resp. Tr) is the initial temperature of the leftmost (resp. rightmost) point x\ (resp. 
Xi) of the domain. A is a parameter that determines the initial stiffness of the temperature gra- 
dient. 

The simulations were performed with the following parameters: the uniform magnetic field 
B?,{t = 0, x\) = 0.001,0.01,0.1,1, the size of the dimensionless domain L = x± — x\ = 5400, 
2 x v max = 12, the ion charge Z = 1, the frequency ratio v e i/uj pe = 0.01, the electron thermal 
velocity such as vth/c = 0.05. The initial electric field is zero over the domain: E\{t = 0, x\) = 
E%(t = 0,xi) = 0. The initial distribution function is a Maxwellian depending on the local tem- 
perature, the density being constant over the domain. The initial temperature profile is chosen 
such as Tl = 0.8, Tr = 1.2 and A = 10. This set of parameters enable us to consider the local 
regime, close to the hydrodynamics (the Knudsen number is about 1/500). The dimensionless 
time step and meshes size are At = 1/500, Axi = L/126, Av = 2v max /32 respectively. The grid 
has 126 points in space and 32 3 points in velocity; 42 processors were used for each simulation 
(CEA-CCRT-platine facility). Domain decomposition on the space domain allows each processor 
to deal only with 3 points in space. The fourth order scheme on the transport part has been used. 

Results are presented in Figures [911111 The typical run time is 24 hours for 40 collision times, 
with that set of parameters. The maximum difference between the numerical and the analytical 
solution are less than 10% for longitudinal macroscopic quantities (heat flux and electric field); 
20% for transverse ones. Transverse quantities have only been considered for simulations presented 
in Figures [10] and [11] where the magnetic field was strong enough so that 

• The establishment of transverse heat flux can be acheived during the simulation time. 

• Transverse quantities cannot be considered negligible compared to longitudinal ones. 

These conditions where fulfilled for B% = 0.1, 1. 
In Figures [9lll[ only results for simulations with B% = 0.001, B3 = 0.1, B3 = 1 are shown, 
respectively. The simulation with B3 = 0.01 proved to show no significance differences with those 
with B 3 = 0.001. 

Results shown Figures [9lfTT1 are revealing an important transient phase before the establish- 
ment of a stationary regime. The oscillations are enforced by the magnetic field, Figure [TT] The 
oscillating electric fields are the consequence of the plasma waves excited by our initial conditions; 
they are damped in a few electron-ion collision times. These plasma oscillations are smoothed out 
by the large time steps we used in simulations, allowed by the implicit treatment of the Maxwell 
equations. However this has a little importance on the asymptotic values and a little impor- 
tance for accuracy. With a larger magnetic field Figure [TTl we observe frequency modulations at 
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Figure 9. Longitudinal (along the temperature gradient) ratios max"^^^^) 

(dashed curve) and j^^ 1 (oscillating curve) are shown against the dimen- 

sionless time. The dimensionless magnetic field is B% = 0.001. Asymptotic be- 
haviour, where the flux is well established, shows good agreement (less than 5% 
error) with analytical solution (Braginskii formalism) , denoted by subscript BR. 



Dimensionless time 

(a) Longitudinal 




Dimensionless time 



(b) Transverse 



Figure 10. Ratios mAXx i^ Fp "> (curve in bold) and maX:c i \^f p \ (dashed curve) are 

max xl (Q B R) 1 7 ma,x xl (E BR ) I / 

shown against the dimensionless time. Longitudinal quantities (along the temper- 
ature gradient) are shown in (a), with about 10% accuracy in the asymptotics. 
Transverse quantities are shown in (b), with about 20% accuracy in the asymp- 
totics. The dimensionless magnetic field is B% = 0.1. 



<^c = Ve,i (corresponding to B% = 1), both on electric fields and heat fluxes. 

In order to investigate Larmor radius effects for simulations presented in Figures [10] and [Til we 
refined the space grid below the dimensionless Larmor radius rj_, = B^ . Therefore, simulation 
presented in FigureHUJhas been done again with the same parameters on the same time period: we 
have refined the grid to 1260 points in space (420 processors). In the same manner, the simulation 
presented in Figure [TT1 has been done again with 6300 grid points in space (2100 processors) and 
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(a) Longitudinal (b) Transverse 

Figure 11. Ratios maX3: i (curve in bold) and maX:c i ^ FP \ (dashed curve) are 

shown against the dimensionless time. Longitudinal quantities (along the tem- 
perature gradient) are shown in (a), with about 5% accuracy in the asymptotics. 
Transverse quantities are shown in (b), with about 20% accuracy in the asymp- 
totics. The dimensionless magnetic field is B% = 1. 



At = 1/1000 (C.F.L. condition), during the same time period. The results prove to be similar 
to those with coarse space grids, both for macroscopic quantities and distribution functions. We 
thus show no dependence on the Larmor radius. Here we remark that the cyclotron frequency is 
always resolved. The time steps are constrained, for most of the cases we treat, by the C.F.L. on 
collision operators. 



6.2.2. Electron transport in the nonlocal regime. The departure of transport coefficients from their 
local values is of interest here. We restrict ourselves to cases where the magnetic field is zero. Then 
it is possible to obtain directly the ratio of effective thermal conductivity to the Spitzer-Harm 
conductivity k/ksh by the relation: 



KSH QSH 

The Spitzer-Harm regime refers to a local regime with no magnetic field. In (|53p . q\ is calculated 
from the numerical solution and qsn from (|48p in the Spitzer-Harm limit. 

Transport coefficients are extracted from the domain where the flux and temperature gradient 
are maximum. 

The wavelength of the temperature perturbation k\ e j in the Fourier space is computed from the 
gradient temperature profile. This enables to obtain a range (due to an uncertainty) for fcA e ,i 
corresponding to this temperature gradient. 

The results will be compared with the analytical formula from |17| 



(54) 
(55) 
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The comparison between the numerical results and the analytical solution are in good agree- 
ment. The three runs have been performed with the same precision for the temperature gradient. 



Parameters 


RUN1 


RUN2 


RUN3 


Size of the domain 


5400 


540 


540 


Stiffness parameter A 


10 


10 


100 


Number of points along the Gradient 


126 


126 


1260 


Number of processors 


42 


42 


420 


Results 


RUN1 


RUN2 


RUN3 




10~ ;i 


0.05 ± 0.03 


0.2 ±0.1 


Analytical k/ksh 


0.998 


[0.93 - 0.67] 


[0.60 - 0.26] 


Numerical k/ksh 


1.03 


0.675 


0.395 



6.3. 2D nonlocal magnetic field generation. We present here results on the nonlocal mag- 
netic field generation during the relaxation of cylindrical laser hot spots, having a periodic repar- 
tition, and for a region of constant density. This stands as a first step to prove the 2D capabilities 
of the solver. The 2D extension of the presented numerical schemes is straightforward on a grid. 

We consider a planar geometry with periodic boundary conditions. For this application, the 
normalizations using collision frequency ()35[) - (|3T[) are used. 

The initial dimensionless temperature profile is T e (x,t = 0) = l + 0.12exp ^— , with R = 5.6. 

We used the following parameters for the simulation: the frequency ratio is v e ^/u pe = 0.003, the 
ion charge Z is assumed to be high, so that we do not consider the electron-electron collision 
operator; here the relaxation only acts with electron-ion collisions on the anisotropic part of the 
electronic distribution function. The electron thermal velocity is such as vth/c = 0.05. These 
parameters are close to those used in [31]. The size of the simulation domain is L = 70 for one 
space direction, 2 x v max = 12 for one velocity direction. Initial electric and magnetic fields 
are zero over the domain. The initial distribution function is a Maxwellian depending on the 
local temperature, the density being constant over the domain. The dimensionless time step and 
meshes size are At = 1/500, Ax = Ay = L/100, Av = 2v max /32, respectively. The grid has 100 2 
points in space and 32 3 points in velocity. 625 processors are used for this simulation. 

The mecanism under consideration here (the magnetic field generation Figure [T2]) . is expained 
in [23], as the results of non parallel gradients of the third and fifth moments of the electronic 
distribution function. We show the magnetic field in Figure [T2l (a) and the cross gradients 

V x ( f fehfdv) X V x ( f fe\H\ 5 dv) in Figure 121(b). 

This mecanism is not due to the magnetic field generation from a Vn e x VT e structure, since the 
density n e remains constant over the domain. 

This structure with eight lobes is the result of the collision operators (of diffusion type) that make 
a particular speckle interact, after a rapid transient phase, with the other surrounding (similar) 
speckles. We note that an important parameter to anayse further such interactions should be the 
size of the speckle over the distance between speckles. 
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Figure 12. Dimensionless magnetic field and cross gradients of high order mo- 
ments (third and fifth) at tv e ^ = 8. 



7. Conclusions 

In the present paper, we have developed high order numerical methods dedicated to plasma 
simulation at a microscopic level. 

A fourth order scheme issued from VFRoe schemes has been introduced in our kinetic context. 
It brings accuracy improvement on the velocity transport term. The second order scheme remains 
interesting for the linear spatial transport term (which is faced to less robustness and accuracy 
constraints) in a 2D, distributed memory context without overlapping between processors (each 
processor communicating with its neighbours only). It involves indeed a reduced stencil allowing 
for a lower minimum number of spatial grid points per processor. 

The Maxwell equations have been discretized with a second order, implicit scheme allowing large 
time steps. We did not find any dependance on the Larmor radius and show that resolving the 
cyclotron frequency is sufficient. The couplings between the equations of the system have intro- 
duced a number of constraints (robustness, accuracy, symmetry) both on the transport scheme 
and the collision operators. Some numerical and physical test cases have validated our approach 
in different regimes of interest for ICF applications, and showed that it is computationally af- 
fordable. We also proposed a validation strategy in the linear regime based on [30], using Green 
kernels. 

Various fundamental studies can be planned on the basis of the actual version of the solver. Col- 
lisional Weibel instability [28], forward and backward collisional Stimulated Brillouin Scattering, 
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studies on the nonlocal interaction beween speckles for plasma-induced smoothing of laser beams 
issues [18], for instance. Also several axis of development are under consideration to bring more 
physics to the model: the ion motion, the extension to regimes relevant to higher laser intensities 
(relativistic regime and large angle collision terms of Boltzmann type) . 
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The relativistic 1D X x 3D V Vlasov-Poisson system extracted from the equations jl])-© reads 



Appendix A. Electrostatic case in the linear regime 
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under the hypothesis: 
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The Vlasov-Poisson can then be set under the following form (transport equation along the space 
directions supplemented by a source term along the v\ direction), after linearization 
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If /W and are periodic and integrable, then their respective normalized Fourier coefficient 
are well-defined. A Fourier series expansion gives Vi > 
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Where L is the size of the domain. The same reconstruction using Fourier series is used for E^ . 
These coefficients verify the following equations, obtained by Fourier transformation performed on 
the equations of the system (j62j) . for all real ki 
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Then introducing the notation f^(t = 0, ki,v) = f( 10 \ki,v), the equation (i64|) can be written 
in the integral form 
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Integrating the equation (|66p over v and injecting in it the relation (|65p . one obtains the following 
integral equation for the density 



(67) 

where 
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These kernels can be computed with the desired accuracy, following [30]. The numerical resolution 
of (|67p finally reduces to the inversion of a triangular linear system. 

Macroscopic quantities such as the density or the heat flux can then be reconstructed using these 
latter equations. 
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Appendix B. Initialisation for the generation of a single X-mode plasma wave 



This test case stands as a validation for the couplings of Vlasov and Maxwell equations. We 
determine initial conditions that trigger a plasma wave at a given wavelength. To do so, Vlasov- 
Maxwell equations are linearized, setting / = /(°) + /, E = E, B = + B around the 
equilibrium state f = , E = 0, B = B^°\ In this appendix, we use the normalization ()23[) - 
([25]) . We assume periodic boundary conditions. The fluctuations of the total pressure tensor are 
neglected with respect to those of the magnetic field. 

Using the conservation law — — h — — = 0, the former hypothesis lead us to solve the system of 

dt OX\ 

six equations with six unknown j\ , j 2 , E\ , E2 , B 3 and n 
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+ 32 , 



Applying time and space Fourier tranform to this system, and identifying Fourier composants 
(n = nexp(— iujt + ikiX\)), the following system is obtain 



-iojji + h+ B^j 2 
-iuj 2 + E 2 - B<®j! 
—icon + ikiji = , 
ikiEi = —h , 
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The dispersion equation of this system reads 
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2 ihB 3 +j 2 , 
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In this equation, the plasma frequency is oj pe = 1 and the cyclotron frequency is tv c = q e m, 
that is ||-B°|| in this dimensionless case. The perturbative term of the distribution function at 
initial time can be determined for a particular solution u of this relation dispersion. 
The Fourier transform is applied on the linearized Vlasov equation 



(72) 
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This equation is expressed in cylindrical coordinates 

v\ = v± cos(V>) , 
v 2 = v± sin(V>) , 

v 3 = v\\ 



where 



Recalling that: 
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with V v fj_ = e±, V v iy, = and V v f|| = eji, where e are vectors in the local basis. Setting 
/(°)(||v|| 2 ) = (2vr)t exp(— ll^f ), and writing 

(vAB).V v /=(V v /Av).B = -B(°)|, with B = (0,0, B^) , 

the kinetic equation ([72]) becomes 

(73) (-iio + cos(^))/ + # (0) f( + / (0) (l|v||>±(£i cos(^) + £ 2 sin(V)) = 0. 

dip 

In order to solve this equation, we decompose the distribution function as a Fourier serie 



imp 



Then from (|73[) 

+00 



2 (-tw + ijfciuj. cos(^) + inB^)f n e int ^ = -/ (0) (||v|| 2 )«j.(^i cos(^) + £ 2 sin(V>)) . 



Multiplying this equation by e im ^, integrating from to 2ir, we obtain 

hoo „ 2 7T 



+ ihv± cos(^) + inB^)f n ^dip 



(74) 



277 



-/ (0) (l|v|| 2 )^± / e^^i cos(V') + E 2 sin(V>))# 

38 



For m = 0, terms are different from zero only for n = — 1, 0, 1. From (|74p comes 



(75) 



hv±f-i - 2w/ + fciui/i = . 



For m = — 1, 



(76) 



ifcivi/o - 2i(w - B(°!)/i + iM±/ 2 = -/o(^>±(£i - i£ 2 ) • 



For m = 1, 

(77) ifciUj_/- 2 - 2i(w + B (0) )/_i + ihvjo = -/(°)(|| V || 2 ) W± (^ 1 + iE 2 ) . 
The case m = — 2 involves /3, 

(78) ifciuj./! - 2(w - 2£?(°))/ 2 + ifci^/a = . 
In the same manner the case m = 2 involves /_3, 

(79) *fci«±/_3 - 2(w + 2fl<°>)/_ 2 + iAsiwj./-! = . 

In order to close the system, the components /_3 and /3 are neglected, and we deduce from 

mm, 



-2(co + 2BW)f„ 2 + ikiv x f- 1 = 0, 

toi/_ 2 -2i(w + #)/_i+iWo = -/(°)(||v|| 2 K(Ei+^ 2 ) , 

kv±f-i - 2a; /o + fci«±/i = , 

ifcWo - 2i(w - B(°))A + ifciu ± / 2 = -/ (0) (l|v|| 2 K(£i - i^a) , 

ifcWi - 2(w - 2B(°))/ 2 = 0. 



The solution of linearized Vlasov equation can be calculted 



( f(t,X,v) = f(°\\\v\\ 2 )+Z+=-oofn(v±)e~ M ^ +m ^ , 

Ex(t,x) = E ie - iujt+iklXl , 

E 2 (t,x) = E 2 e- iult+iklXl , 

, B(t,x) = B(°) + B 3 e~ iu;t+iklXl . 



The dispersion relation (|7T1) provides with a particular Then we obtain the following results 
for the construction of the initial solution, 



/(0,x, W ) = /(°)(||v|| 2 )+ U(v±)e 



ikix+inip 



n=-2 



With the expressions 

= i(-Aoo z Ei - 4 iuo 3 E 2 + 12 iuj 2 B^ E 2 + 12 lo 2 B^ E Y -SWB^fuEh 

+ fci V 2 w + 3 ik! 2 v ± 2 uj E 2 - 8i\\B^ feu E 2 - 4 ik 1 2 v± 2 B^ E 2 )v± 2 h, 

= 2 iv ± (Ex k^v^oj 2 + 4 iB^ uj 3 E 2 - 16 ||£(°> fu E x - 16 fw £ 2 

+ 3iE 2 k^v^uj 2 — AEiuj 4 — 8iE 2 fci V 2 ||B (0) || 2 + 2k 1 2 v ± 2 B^ lo E x 

+ 2 ik! 2 v ± 2 B^ u E 2 + 16 ^1 || V + 16 i£ 2 ||B(°) feu 2 + 4 B< > u 3 ^ 

- 4iE 2 u 4 ), 

= 2 iv± 2 h{16 ||S(°) + fci V 2 w ^1 - 4w 3 ia + 4iu 2 B^ E 2 

- 16i||B(°) H 3 ^ + 2ifciV 2 £(°> £2), 

= 2 i(-2 + (-4 ik! 2 v ± 2 BW e 2 + fci V 2 w A - 3 ik 1 2 v± 2 u E 2 

- 12uj 2 B^ E x + \2iuj 2 B^ E 2 -Aoj 3 Ei+Aiuj 3 E 2 -S\\B^\\^ujEi 
+ 8z||^W|| 2 c^ 2 ), 

= ik 1 v 1 _ 2 (-4iki 2 v± 2 B^ E 2 + k^v^uj E! - 3ik! 2 v± 2 uE 2 - 12u 2 B^ E x 

+ l2ico 2 B^ E 2 -4lu 3 E 1 + 4iu 3 E 2 - 8\\B^fuj E t + 8i\\B^\\ 2 u E 2 ), 
where 

D = uj (64 ||B^ || 4 - 16 fci V V + W + 16 fci V 2 ||B (0) f + 3 fci V 4 - 80 ||B(°> || V). 

k\v±_ being small with respect to B^ and cj, powers of k\v±_ can be neglected compared to these 
terms. The solution can be written 

,{~ 2 |ox A = w ± 2 A;i(-4cj 3 ^i -4^ 3 J B 2 + 12icj 2 J B(°)^2 

+ 12a; 2 B( ) A-8||5(°)|| 2 wSi - 8i\\B^\\ 2 u E 2 ), 
= 2iv ± (4iB(°) w 3 £ 2 - 16 ||£(°)|| 3 u;£i - 16i||B(°) £ 2 -4£iw 4 



/ (0) (l|v| 


\ 2 )D 


7-1 




/ (0) (l|v| 


\ 2 )D 


fn 




/ (0) (l|v| 


\ 2 )D 


fl 




/ (0) (l|v| 


\ 2 )D 


A 




/ (0) (l|v| 


\ 2 )D 



/ (0) (l|v| 


\ 2 )D 


f-i 




/ (0) (l|v| 


\ 2 )D 


/o 




/ (0) (l|v| 


\ 2 )D 


fl 




/ (0) (l|v| 


\ 2 )D 


A 




/ (0) (l|v| 


\ 2 )D 



+ 16 E x II V + 16 iE 2 ||S(°) ||V + 4 w 3 £i - 4 w 4 ), 

= 2w ± 2 A;i(16|| j B( )|| 2 lj j Bi + k 1 2 v 1 _ 2 ioE 1 -4w 3 li + 4iuj 2 B^ E 2 

- leipwn 3 ^), 

= 2iv± (u -2BW)(-12lu 2 B ( °1 E x + 12iu 2 B^ E 2 - 4<J 3 E 1 + 4iuj 3 E 2 

- 8\\B^\^E 1 +8i\\B^\fujE 2 ), 

= ik lV± 2 (-12 w 2 B(°» £1 + 12kj 2 5(°) E 2 - 4w 3 £i + 4iu; 3 £ 2 

- 8||S(°)|| 2 L l j J Ei + 8i|| J B( )|| 2 cj J E2), 

40 



where 



D = u (64 \\B^ || 4 + 16a; 4 - 80 



,2 



.2\ 



We choose to initialise the perturbation from the amplitude of the magnetic field: 

B 3 = A where A £ [0, 1]. 

Then from the system ()70|) and the dispersion relation (|7ip . we deduce the values of E±, E2 and 
thus reconstruct the fi, 



iB 3 ( u 4 (3 2 - uj 2 h 2 - uo 2 (3 2 - || V/3 2 + ||B(°) || V 



-^1 — ; — ^ „/ n N , -E 1 : 



A;i/3 2 B(°) ' z fci 
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